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We study the Bose-Einstein condensation phase transition in a weakly in- 
teracting gas through a perturbative analysis of finite systems. In both the 
grand canonical and the canonical ensembles, perturbation theory suffers from 
infrared divergences and cannot directly determine the transition temperature 
in the thermodynamic limit. However, in conjunction with finite size scaling, 
perturbation theory provides a powerful calculation tool. We implement it 
here to estimate a shift in the transition temperature in the canonical ensem- 
ble consistent with grand canonical calculations. 
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Wilkens, Illuminati, and Kramer [1] derive a surprising and extremely 
interesting result, namely that as a function of interaction strength the 
transition temperature of a dilute Bose gas behaves differently in the canon- 
ical and grand canonical ensembles. They conclude that in the canonical 
ensemble, the transition temperature decreases with increasing interaction 
strength, whereas the opposite behavior has been established in calcula- 
tions in the grand canonical ensemble [2, 3]. The latter agree with nu- 
merical calculations in the canonical ensemble which extrapolate to the 
thermodynamic limit using finite size scaling techniques [4, 5] . 

In this paper we reformulate the ideas of Wilkens et al. [1] to provide a 
more explicit comparison with its grand canonical counterparts and with 
numerical calculations. The calculation is perturbative in nature; we show 
that such a perturbative scheme fails in the thermodynamic limit due to the 
presence of infrared divergences. As in the grand canonical ensemble, long 
distance properties dominate thermodynamic quantities of the system, as 
expected for a second order phase transition. Standard finite size scaling, 
as used in [5], provides a tool to overcome these difficulties, giving a precise 
scheme to extrapolate to the thermodynamic limit. 
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This paper is structured as follows. First, to lowest order in the in- 
teraction strength, we derive a perturbative expression for the transition 
temperature. Next, by considering higher order terms, we show that the 
perturbation series diverges in the thermodynamic limit. Finally, we use 
finite size scaling in conjunction with the perturbation series to correctly 
derive the shift in the transition temperature. Our first section closely fol- 
lows [1], as one of our primary goals is to reconcile their calculation with 
conventional results. 

Consider a system of bosons interacting through a short range potential 
described by the Hamiltonian 



q 

where 



H int = ^^J2 b t b %- kb P+ k ' ( 2 ) 

pkq 

and a is the scattering length, b q is the annihilation operator for particles 
with momentum q, m is the particle mass, and V = L 3 is the volume of 
the system. The free energy to first order in a is 

F(a)=F Q + (H int ), (3) 

where the expectation value is in the free ensemble, and F = F(a = 0) is 
the free energy of the non-interacting system. The basic technical difference 
between the canonical and grand canonical ensembles arises in factoring the 
four point expectation value in Eq. (3); in the grand canonical ensemble 

E (<*&-*>$W + (blb p+k )(blb q _ k )) (4) 

pqk pqk 

= 2<A0 2 , (5) 

where (N) = ^2 k {b k bk) is the mean number of particles. In the canonical 
ensemble 

EWVfcWc = 2N " Y,( N p( N p + !))' ( 6 ) 

pqk p 

where N p = b^bp. Thus the first order shift in free energy differs in the two 
ensembles by the term 

(H int ) c ~ <tfi„t>GC = -^£<Wp + 1)>- (7) 

p 
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Although the coefficient of this sum is of order 1/V, the sum itself diverges 
as V 4 / 3 at the Bose condensation transition. This term of order V 1 / 3 is 
the source of the shift in the critical temperature described in [1]. The 
contribution proportional to ./V 2 can not affect the critical temperature, 
since this term just adds a constant to the free energy. We now restrict 
ourselves to the canonical ensemble to determine the role of the term (7). 
In all subsequent manipulations both the number of particles iV and the 
system size V are held fixed. 

The order parameter for the Bose-Einstein condensation transition is 
No, the number of condensed particles. To explore the transition we study 
the probability of having N particles in the condensate. This distribution 
function, P(N ), is 

P{Nq) = _L Tr e^ = ie-W»), (8) 

Zjjq N,N fixed Zjv 

where the trace is taken at fixed N and No, and (3 is the inverse tem- 
perature. This equation defines the free energy F(N,N ), and uses the 
partition function, Zn — J2n e~ /3F ^ N ' N °\ to normalize the probability 
distribution. The high temperature normal phase of the system is charac- 
terized by a monotonically decreasing P(Nq), while the low temperature 
condensed phase has P(Nq) peaked at iVo ^ 0. At sufficiently high temper- 
atures, No is Gibbs distributed, P(N ) oc e fi ^ Na , where fj, = dF/dN < 
is the chemical potential. At zero temperature in the non-interacting gas 
all the particles are condensed and P(N ) = Sn 0j n- At some intermediate 
temperature the distribution becomes flat at No — 0. Wilkcns et al. de- 
fine T c by extrapolating this crossover temperature to the thermodynamic 
limit. 

In terms of F(N,No), Wilkens et al.'s criterion for T c gives an implicit 
equation for the critical temperature of the interacting system Ti°\ 



dF(N,No) 



dN 



Ar„=0, T=T c (o) 



= 0. (9) 



As in Fermi liquid theory, 8F/dN is the energy of a k — quasiparticle 
measured from the chemical potential [6], and can therefore be expressed 
as 

OF 

_ = S(fc = 0, W = 0)- M , (10) 

where Yi{k,uj) is the self-energy at momentum fc and energy lo. Thus this 
criterion for the critical temperature is essentially that used by Baym et 
al. [3] in the grand canonical ensemble. An important difference between 
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the two approaches is that in the present calculation only quantities at 
iV = arc involved. In the canonical ensemble the fluctuations in No are 
very large at the critical temperature, (Nq) — (No) 2 ~ N 4 / 3 [7]; as we shall 
see, the criterion (9) yields a qualitatively different shift in the transition 
temperature if the derivative is evaluated at the expectation value of No 
rather than at No = 0. 

We now expand the criterion (9) for T c in powers of a, to calculate 
perturbatively the transition temperature, , of the interacting system. 
Since dF/dNo is evaluated at T — Tc a \ we must consider not only the 
explicit variation of F with a, but also the implicit contribution due to the 
dependence of T on a. We use the decomposition F(a) = Fo + AF(a), 
where Fo is the free energy of the non-interacting gas and AF(a) is the 
correction due to interactions. In the free system, the condensate only 
contributes to the free energy by reducing the occupation of other modes, 
i.e., 

dF (N,N ) dFo(N,N ) 

dN ~ ON = ~ M0 ' (U) 

which defines the free chemical potential a function of N, N and T. 
By construction, when No — this chemical potential vanishes at the 
transition temperature of the non- interacting gas, Tc°\ and to first order 
in the interaction, 



dF 



8N 



jV =0 : T=T e <a) 



M o(^o = 0,T = T c (°)) (12) 
d[io 



-AT C 

dT 



]V o =0,T=T c (o) 



O(a'). (13) 



The derivative is taken at fixed N and N Q , and AT C = T c (a) - T c (0) is the 
shift in the transition temperature for scattering length a. Thus, to first 
order in a, the left hand side of Eq. (9) becomes d(AF)/dN — AT c d^o/dT, 
evaluated at T = T c (0) and N = 0. Solving for AT Cl we have 



AT _ (d(AF)/dN Q ) 



(d^o/dT) 



N =0,T=T^ 



_rp(0) 



(14) 



Aside from the use of continuous derivatives in place of Wilkens et al.'s 
discrete derivatives, this is the result of [1]. Correctly evaluating these 
functions for a finite sized system is challenging. We estimate their magni- 
tude by replacing the canonical expectation values in Eq. (6) by the grand 
canonical result (Nk(Nk + l)) = 2(Nk}((Nk} + l), and approximately writing 
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(Nk) ~ (e /3 ^ fc_AiQ ' — This assumption provides a simple relationship 

between N = N a + ^2 k {N k ) and u - Introducing the number of excited 
particles A cx — N — N , we may write 



d(AF) d(AF) _ fd(AF)\ ( <9 Mo 



8N 9iV cx V 9^0 J \dN 



( d(AF) \ ( g Mo 
V d^o J \dN cx 



(15) 



Since all quantities are evaluated in the free ensemble, the derivatives are 
straightforwardly evaluated, leading to 

AT C 8irh 2 a ^ 3 

where (Nk) is evaluated at /l*o = 0. The sum is infrared divergent, scal- 
ing as V 2 , and yielding a finite temperature shift proportional to —an 1 / 3 , 
where n = N/V. The constant of proportionality is of the same order 
of magnitude as the one calculated in [1] using a sophisticated series of 
asymptotic expansions; its exact numerical value is unimportant here. The 
key observation is that contrary to the expected behavior, the temperature 
shift predicted by this argument is negative. 

This negative temperature shift depends crucially upon the constraint 
N a = 0. At finite N , the numerator of Eq. (14) has an additional contri- 
bution due to the derivative of the p = term of the sum in Eq. (7). This 
contribution has the opposite sign, and dominates when Nq <~ iV 2 / 3 , yield- 
ing a positive temperature shift. As already emphasized, at the critical 
temperature, the expectation value of No is of order iV 2 / 3 . 

We now explore the validity of this perturbation expansion, demonstrat- 
ing that it breaks down in the thermodynamic limit. Higher order terms 
in the expansion of the free energy (3) involve higher powers of the in- 
teraction Hi n t. As in first order, the most divergent terms occur when 
all of the momenta are equal and, at T c , these terms are of relative size 
((/3H int ) m ) J ' {(l3H int ) m - 1 ) ~ aL/X 2 , where L = V 1 ^ is the length of the 
system, and A = {2-kTi 2 jmk^T) 1 ! 2 <~ n~ 1 l z is the thermal wavelength, of or- 
der the interparticle spacing. Thus we see that the perturbation expansion 
is valid only for sufficiently small i] = aL/X 2 . 

In any finite system the transition temperature calculated above corre- 
sponds to a crossover of finite width ST. One can estimate the width of 
the crossover from the fluctuations in the number of condensed particles; 
since the latter scale as [7], 

ST/T ~ SN /N ~ A~ 1/3 . (18) 
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As long as ST < AT C , the shift is well defined. Since AT oc -an 1 / 3 , the 
ratio AT /ST is of order aL/X 2 = 77, and the shift only is well defined for 
r\ 3> 1. Thus in the limit of small 77, where the expansion of the free energy 
converges, the calculated change in T c is smaller than the width of the 
transition and cannot be physically significant. In the other limit i)> 1 
the expansion of the free energy breaks down. Hence this calculation, as 
it stands, cannot tell us anything about the transition temperature of a 
weakly interacting Bose gas. 

Although perturbation theory breaks down in the the thermodynamic 
limit, it can be used to learn the properties of small systems where L <C 
A 2 1 a. We now discuss how finite size scaling enables one to learn about the 
L = 00 phase transition by investigating how physical quantities scale with 
L in these small systems. This technique is commonly used in numerical 
simulations where it is not feasible to simulate an infinite system. The 
central assumption of finite size scaling is that sufficiently close to the 
critical point all physical quantities scale as functions of the ratio of the 
correlation length £ to the system size L. For example the order parameter 
scales as 

M „ L-»$(L/0, (19) 

where y — (5 jv = 1 is the ratio of the critical exponents for Nq/V and the 
correlation length, and <J> is a scaling function. As L/£ — > 00, this function 
must diverge as 4>(L/£) <~ {L/£) y , while as £/£ — > 0, $ approaches a 
constant. The latter limit gives a systematic method for finding the critical 
point (£ — » 00) by looking solely at the properties of a finite system. In 
numerical calculations [4, 5] one plots L V (N )/V, or a related quantity 
such as the superfluid density, as a function of temperature for different 
system sizes. According to the scaling hypothesis all of these curves should 
intersect at the critical temperature. 

Relations similar to Eq. (19) also hold for higher moments of the order 
parameter and imply that at T c the probability distribution function can 
be written as 

P(N ) = ^(N X 2 /L 2 ), (20) 

with some scaling function ip. We calculate the critical temperature of the 
interacting system by finding the temperature at which P(Nq) has this 
scaling form. 

Before using this procedure to calculate AT C , we verify that the scaling 
relations, Eqs. (20) and (19), hold in the non-interacting gas. We first derive 
Eq. (19) in the grand canonical ensemble where the argument is particularly 
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simple. The general strategy is to fix the average density n = N/V and the 
temperature T, and look at how the order parameter iV varies with the 
system size L. To carry out this approach we need an expression for the 
chemical potential /i as a function of n, T, and L, which requires inverting 
the relationship, 



V eP^u-ii) _ i ( 21 ) 

k 

« ^ff3/2(e^) + JVo/V, (22) 

where 33/2(2) = ^ . z J /j 3 / 2 is a polylogarithm function. The inversion is 
performed by expanding (22) in powers of noting that No ~ —l/(3fj,. 
One finds 



1 



C(3/2 ) + A(-^_2^^) + . 



; (23) 



for <~ X 2 /L 2 <~ 1, the neglected terms are of relative order A/L. The 
terms proportional to and \J—$\i are respectively the contributions 

from the condensed and non-condensed particles. Finding as a function 
of n, T, and L, requires solving a cubic equation. We define the function 
F(x), plotted in Fig. 1, as the solution to 

^—-2^F(x)-x = Q, (24) 

so that the chemical potential can be expressed as, 

fa = -^F^(A 3 n-C(3/2))). (25) 

The positive, monotonic F(x) has the properties 

F(0) = (Air)- 1 / 3 (26) 
F(x) — > x 2 /4tt (27) 

X— ► — OO 

F(x) — > l/x. (28) 

x— >+oo 



Thus, as L — > 00, the order parameter iV ~ —l/(3fj,, has three distinct be- 
haviors, corresponding to non-condensed, critical, and condensed regimes, 
depending on whether n is less than, equal to, or greater than C(3/2)/A 3 . 
In the non-condensed regime, N is microscopic, in the condensed regime, 
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F(x) 




-5 5 10 15 20 X 

FIG. 1. The function F(x), defined as the solution to Eq. (24), which relates the 
chemical potential and density of a non-interacting ideal bose gas via Eq. (25). 

N is extensive, and at the critical point, N scales as Eq. (19), with the 
predicted exponent, y = 1, i.e., 

n < C(3/2)/A 3 , N ~ L° (29) 
n = C(3/2)/A 3 , iV ~ L 2 (30) 
n > C(3/2)/A 3 , N Q ~ L 3 . (31) 

The scaling at T c , Eq. (30), is consistent with (18), since at the critical 
point the mean value of No is of the same order as the fluctuations 6 No. 

We now perform the equivalent calculation in the canonical ensemble. 
Starting from the expression N = N +J2k(Nk), and approximating (Nk) — 
^ e 0(,e k -no(N o )) _ we solve for no as a function of No- Integrating 

Ho{N ) with respect to N yields F(N,N ), and P(N ) cx e - 0F( - N ' N °\ As 
P(No) is peaked about the point where /i = 0, we can expand J2k(Nk) in 
powers of ^o- Standard asymptotics [8] yields a series which converges for 
|/?Mo| < n\ 2 /L 2 , 

N CK = ^(N k )=^am+[j)\[^)+0(L/X),(32) 

OC fa 

f& = £tT C 3(* + 1), (33) 
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TABLE 1. 

Values of the lattice sum Cs(k), defined in Eq. (34). 



k=l 2 


3 


4 


5 


6 


-2.8374 1.6752 


0.5419 


0.4278 


0.5039 


0.7741 



where the coefficients Cd(fc) are sums of the form, 

n n u t?...,n d (n 2 1 +n 2 + --- + n 2 ) 

the prime denotes that the term m = n 2 = • • • = rid = is omitted. These 
constants are tabulated in Table 1 . Inversion of the series gives 

f-TT^M - -g#^ 2 + • • ), (35) 



L 2 V C 3 (2) 2(C 3 (2)) 2 - 
M = iVo^J + C 3 (l) - ^ - ^C(3/2)) , (36) 

from which the free energy is 

f3F (N,N ) = - JdN ((3rio) (37) 

^« + ^) M2 + » M3 + ''' (38) 
= [3F (N)+g(M). (39) 

where Fo(N) is an extensive function which is independent of No, and g(M) 
is defined by Eqs. (38,39). If n = N/V = C(3/2)/A 3 , then M depends on 
L only through the variable N$\ 2 /L 2 , implying that P(Nq) is of the form 
(20). For any other value of the density, M has additional L dependence, 
and P(N ) does not have the desired form, implying that the scaling (20) 
holds only at T c . In Fig. 2 we plot P(Nq) at the critical point. The 
Gaussian approximation, where only the term proportional to M 2 is kept, 
is also plotted, and agrees quite well with the full result. 

Having verified the scaling relationships for P(Nq) in the non-interacting 
gas, we now search for the critical temperature of the interacting gas by 
perturbatively calculating P(Nq) and finding the temperature at which 
scaling holds. We continue to use the approximations that {Nk(Nk + 1)) = 
2(N k )((N k ) + 1), and (N k ) w ( e ^ e <=-«>) - l)" 1 . To first order in a, at 
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1 2 3 4 5 6 7 L 2 

FIG. 2. Probability distribution P(Nq) for the number of condensed particles at 
the critical temperature in a non-interacting Bose gas within the canonical ensemble. 
Dashed line is a Gaussian approximation. At the critical point, P has the scaling form 
(20). 



temperature T = + AT, the free energy is, 



AT A 2 

0F[N, N ] = PF + g(M ) - — j^N g'(M ) (40) 



+^I c(W(Mo) " St 1 ( *° + 2 E + !) 



where the argument M = N (X 2 /L 2 ) + C 3 (l) is the scaled condensate 
number measured from the peak of the distribution. The first two terms are 
the free energy of the non-interacting gas at Tj ^ , while the remaining terms 
give the first order corrections in AT and a. These corrections are only 
small if f] = aL/\ 2 <C 1. The sum J^k^o ^k(Nk + 1) can be identified with 
cW ex /<9(/3^o), and can be expressed as a series in P/j, L 2 /X 2 via Eq. (32). 
Using Eq. (35) to eliminate fi, the corrections are 

AT A 2 

F(T) = F (T C <°>)-— ^JW(Mo) (41) 



L 

+ X 



3C(3/2) AT (I jr C 3 (3) „, 2 



T V^ 3 (2) 2C 3 (2) 
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+ (l + £M2^W ) M , + ...)]. 

Comparing with Eq. (20), we see that scaling holds if and only if the factor 
multiplied by L/X vanishes. Eliminating the coefficient of the first power 
of M enforces scaling near the peak of P(N n ), in which case 



-T= - 3^(3/2) ^(1)^3(2) + C 3 (3)) (42) 
w 1.6em 1/3 . (43) 

The coefficient 1.6 should be compared with the numerical value 2.3 cal- 
culated by Holzmann and Krauth [5]. The discrepancy lies within the 
accuracy expected of our approximations. The important point to note is 
that the coefficient is positive and of order unity. 

The neglect of terms of higher order in a during the calculation is based 
on the assumption that they do not change the structure of the scaling 
function. (We note that recent calculations of (j> 4 theory on a lattice [9, 10], 
may indicate that this assumption is not valid.) A more involved study, 
where these higher order terms are explicitly calculated would help verify 
whether perturbation theory is valid within finite size scaling. 

At this point it would be appealing to repeat the above calculation in 
the grand canonical ensemble and explicitly verify that the two ensembles 
yield the same shift in the transition temperature. In the grand canonical 
ensemble, first order perturbation theory changes the energy of each mo- 
mentum state by the same amount. This shift can therefore be absorbed 
into the chemical potential, leaving the transition temperature unchanged. 
The first effects start at higher order; exploring how higher order perturba- 
tion theory in conjunction with finite size scaling can be used to calculate 
the shift of T c in the grand canonical ensemble will be discussed in a future 
publication. 

In summary, we demonstrate that infrared divergences prevent the direct 
application of perturbation theory to calculating the transition temperature 
of a dilute Bosc gas in the canonical ensemble. We use scaling arguments 
to circumvent this problem and to evaluate AT C within the canonical en- 
semble, finding results which are consistent with numerical calculations, 
and with analytic results based on the grand canonical ensemble. 



12 



EJ MUELLER, G BAYM, AND M HOLZMANN 



ACKNOWLEDGMENT 

This work was supported in part by NSF Grant PHY98-00978 and the NASA Micro- 
gravity Research Division, Fundamental Physics Program. We are grateful to Fabrizio 
Illuminati and Werner Krauth for seminal discussions. 

REFERENCES 

1. Wilkens M, Illuminati F, and Kramer M 2000 J. Phys. B: At. Mol. Opt. Phys. 33 
(2000) L779. 

2. Baym G, Blaizot J-P, Holzmann M, Laloe' F and Vauthcrin D 1999 Phys. Rev. Lett. 
83 1703. 

3. Baym G, Blaizot J-P, Zinn-Justin J 2000 Europhys. Lett. 49, 150. 

4. Griiter P, Ceperley DM, and Laloe F 1997 Phys. Rev. Lett. 79, 3549. 

5. Holzmann M and Krauth W 1999 Phys. Rev. Lett. 83, 2687. 

6. Baym G and Pcthick C, 1991 Landau Fermi liquid theory: concepts and applications, 
J Wiley and Sons, New York. 

7. Illuminati F, Navez P, and Wilkens M 1999 J. Phys. B: At. Mol. Opt. Phys. 32, 
L461. 

8. Ziflf R, Uhlenbcck G, and Kac M 1977 Phys. Rep. 32, 169. 

9. Arnold P and Moore G 2001 cond-mat/0103228. 

10. Kashurnikov V, Prokof'cv N, and Svistunov B 2001 cond-mat/0103149. 



